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Abstract 

Improved algorithms for the solution of the three- 
dimensional time-dependent Euler equations are presented for 
aerodynamic analysis involving unstructured dynamic meshes. 
The improvements have been developed recently to the spatial 
and temporal discretizations used by unstructured grid flow 
solvers. The spatial discretization involves a flux-split 
approach which is naturally dissipative and captures shock 
waves sharply with at most one grid point within the shock 
structure. The temporal discretization involves either an 
explicit time-integration scheme using a multi-stage Runge- 
Kutta procedure or an implicit time-integration scheme using 
a Gauss-Seidel relaxation procedure which is computationally 
efficient for either steady or unsteady flow problems. With 
the implicit Gauss-Seidel procedure, very large time steps 
may be used for rapid convergence to steady state, and the step 
size for unsteady cases may be selected for temporal accuracy 
rather than for numerical stability. Steady flow results are 
presented for both the NACA 0012 airfoil and the ONERA M6 
wing to demonstrate applications of the new Euler solvers. The 
paper presents a description of the Euler solvers along with 
results and comparisons which assess the capability. 


Introducti on 

Considerable progress has been made over the past two 
decades on developing computational fluid dynamics (CFD) 
methods for aerodynamic analysis. 1 ^ Recent work in CFD has 
focused primarily on developing algorithms for the solution of 
the Euler and Navier-Stokes equations. These methods of 
solution typically assume that the computational grid has an 
underlying geometrical structure. As an alternative, 
algorithms have been developed recently which make use of 
unstructured grids . 8 * 1 6 j n two dimensions these grids are 
typically made up of triangles and in three dimensions they 
consist of an assemblage of tetrahedra. The unstructured grid 
methods have been demonstrated to have distinct advantages 
over structured grid methods in that they can easily treat the 
most complex of geometric configurations and they also easily 
allow for adaptive mesh refinement to treat very complicated 
flow physics. Both of these advantages result from the 
flexibility of orienting and numbering the elements which 
make up the mesh in an arbitrary fashion. Many of the 
methods based on unstructured grids, however, use a spatial 
discretization based on central differencing with explicit 
artificial dissipation, and use temporal discretizations 
involving explicit time-marching such as a multi-stage 
Runge-Kutta time integration. The explicit artificial 
dissipation used in such schemes tends to smear shock waves 
over several grid cells and requires the tuning of free 
parameters that scale the dissipation. Also, the explicit time- 
integration has a step size that is limited by the Courant- 
Friedricks-Lewy (CFL) condition to very small values. 
Consequently, thousands (and occasionally tens of thousands) 
of time steps are required to obtain steady-state solutions, and 
thousands of steps per cycle of motion are required for 
unsteady solutions. Therefore, the purpose of the paper is to 
report on improvements that have been developed recently to 
the spatial and temporal discretizations of the unstructured 
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grid flow solvers which resolve the numerical issues 
described above. The spatial discretization now involves a so- 
called flux-split approach based on the flux-vector splitting 
of van Leer. 17 The flux-split discretization accounts for the 
local wave-propagation characteristics of the flow and 
captures shock waves sharply with at most one grid point 
within the shock structure. A further advantage is that the 
discretization is naturally dissipative and consequently does 
not require additional artificial dissipation terms or the 
adjustment of free parameters to control the dissipation. 
Furthermore, in addition to an explicit time-integration 
scheme involving a multi-stage Runge-Kutta procedure, the 
temporal discretization can alternatively be based on an 
implicit time-integration scheme involving a Gauss-Seidel 
relaxation procedure. This relaxation procedure is stable for 
very large time steps and thus allows the selection of the step 
size based on the temporal accuracy dictated by the problem 
being considered, rather than on the numerical stability of the 
algorithm. Consequently, very large time steps may be used 
for rapid convergence to steady state, and an appropriate step 
size may be selected for unsteady cases, independent of 
numerical stability issues. The work reported herein 
represents the extension to three dimensions of some of the 
methods developed originally by the author in two 
dimensions. 18 In Ref. 18, steady and unsteady flow results 
demonstrated the accuracy and efficiency of the improvements 
to the spatial and temporal discretizations for the NACA 0012 
airfoil. To assess the performance of the three-dimensional 
algorithm, calculations were performed for the ONERA M6 
wing 19 at a freestream Mach number of M^ - 0.84 and an 
angle of attack of a - 3.06*. This case is an AGARD standard 
case for the assessment of inviscid flowfield methods, 19 where 
experimental steady pressure data are available for 
comparison with calculated pressures. 


Eutec-Equat ion s 

In the present study the flow is assumed to be governed by 
the three-dimensional time-dependent Euler equations which 
may be written in integral form as 

^ Jn QdV + J 3n (En « + Fr V + Gn *> dA = 0 (1) 

where Q is the vector of conserved variables representing 
mass, momenta, and energy, and E, F, and G are the convective 
fluxes in the x, y, and z directions, respectively. The second 
term is a boundary integral resulting from application of the 
divergence theorem. 


Spatial Discretization 

The spatial discretization is based on Van Leer's 17 flux- 
vector splitting which is herein implemented as a cell- 
centered scheme where the flow variables are stored at the 
centroid of each tetrahedra and the control volume is simply 
the tetrahedra itself. Consequently the spatial discretization 
involves a flux balance where the fluxes along the four faces of 
a given tetrahedra are summed as 

£ HA » =£ T(EA x +FA y +GA z ) (2) 
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where T is a transformation matrix which is required to rotate 
the fluxes into a locally Cartesian coordinate system that has a 
coordinate direction normal to the face. In Eq. (2) A x , Ay, and 
A z are the directed areas of the face in the x, y, and z 

coordinate directions, respectively, and A*=A* + Ay+A*. 

In general, the flux vector H is split in a one-dimensional 
fashion into forward (H+) and backward (H-) vectors as 

H - H + (q‘) + H-(q+) (3) 


similar to Eq. (6) is also used. Futhermore, in calculations 
involving upwind-biased schemes, oscillations in the solution 
near shock waves are expected to occur. To eliminate these 
oscillations flux limiting is usually required. The flux limiter 
modifies the upwind-biased interpolations for q* and q + such 
that, for example 

<T = Qj +7l(1-KS)A_ +(1 + KS)A + ] (8) 
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where the notation H + (q _ ) and H’(q + ) indicates that the 
fluxes H ± are evaluated using upwind-biased interpolations of 
the primitive variables q. In two-dimensions for a given 
triangle j, for example, and considering the diagram in Fig. 
1{a), the upwind-biased interpolation for q* along the edge 
between triangles j and k is defined by 


where s is the flux limiter. In the present' study, a 
continuously differentiable flux limiter was employed which 
is defined by 


2 a. A+ + e 

A?. + A* + £ 


(9) 


q" = qj + -[(l-K)A_+(1+K)A + ] 


where e is a very small number to prevent division by zero in 
( 4 ) smooth regions of the flow. 


where A + = q k - qj 

a- = q, - 


(5a) 

(5b) 


In Eqs. (4) and (5), q i and q k are the vectors of primitive 
variables at the centroids of triangles j and k, respectively, 
and q^ the vector of primitive variables at node i, is 
determined by an average of the flow variables in the triangles 
surrounding node i. The upwind-biased interpolation for q + 
along this edge is determined similarly using the flow 
variables at centroids j and k and the flow variables at node I. 
In three-dimensions a similar interpolation procedure is used 

where q^ and q k are the vectors of primitive variables at the 

centroids of tetrahedra j and k, respectively, and q ( , the vector 
of primitive variables at node i, is determined by an average of 
the flow variables in the tetrahedra surrounding node i. 

The parameter k in Eq. (4) controls a family of difference 
schemes by appropriately weighting A_ and A + . On structured 
meshes it is easy to show that k = -1 yields a fully upwind 
scheme, k - 0 yields Fromm's scheme, and k - 1 yields central 
differencing. The value k = 1/3 leads to a third-order- 
accurate upwind-biased scheme, although third-order 
accuracy is strictly correct only for one-dimensional 
calculations. Nevertheless, k - 1/3 was used in the calculations 
presented herein. 

On highly stretched meshes in two dimensions, the formula 
for A + is modified to be 


A- = (Qk ~ Qj) 

a + b 1 


(6) 


where a and b are the distances from the midpoint of the edge to 
the centroids of triangles j and k, respectively, as shown in 
Fig. 1(b). This formula weights the flow variables in the 
interpolation formula (Eq. (4)) differently to account for the 
stretching of the mesh. For example, by substituting Eq. (6) 
into Eq. (4) and letting k - 1 yields 


b a 

q = — -qj + — -qk 

a + b a + b 


(7) 


For the case shown in Fig. 1(b), Eq. (7) clearly gives more 
weight in the calculation of q* to the flow variables at centroid 
j than to the flow variables at centroid k, since b>a. On highly 
stretched meshes in three dimensions, a modified formula 


• centroid 



( a ) centroids and nodes used in construction 
of upwind-biased flow variables. 



( b ) distances between centroids and midpoint 
of edge used in Eqs. (6) and (7). 
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Fig. 1 Diagrams illustrating details of the flux-split Euler 
algorithm implementation. 



Temporal Discretization 


The unsteady Euler equations may be Integrated in time 
using either explicit or implicit time-marching procedures. 
Each of these procedures is described briefly In the following 
subsections. 

Explicit Ti me- M arching 


X T " 1 l H+ (q')+H-(q + )r 1 A, > +H "(<J + >lX 

-KXT-VA.JAQj + f-r’A-A.AQn, 

m-1 

(14) 


The Euler equations are integrated in time by assuming 
that the conserved variables represented by Q are constant 
within a control volume which yields 

~ (V Q ) + C (Q) = 0 (10) 

dt 

where C represents the convective fluxes and V is the volume 
of a given tetrahedron. These equations are integrated in time 
using an explicit four-stage Runge-Kutta time-stepping 
scheme given by 

Q (0) = Q n 

Q< 1 > = Q(0)_i AIc/qW, 

4 v 

q<2> = qW.^cjqO)) 

Q<3) = Q(0)_^£t C(Q (2), (11) 

q<4) = q( 0) _ At C{Q (3>, 

V 

Q" + 1 = Q< 4 > 

To accelerate convergence to steady state, implicit residual 
smoothing and local time-stepping are used. The residual 
smoothing allows the use of CFL numbers that are larger than 
that dictated by the stability of the original scheme. This is 
accomplished by averaging implicitly the residual with values 
at neighboring grid points. These implicit equations are solved 
approximately by using several Jacobi iterations. The local 
time-stepping uses a maximum allowable step size at each grid 
point as determined by a local stability analysis. 

Implicit Time-Marching 


In this equation the last summation on the right hand side 
involves AQ m , the change in the flow variables in the four 
tetrahedra adjacent to tetrahedron j. Also, the exact jacobians 

A + and A* are determined by differentiation of H* and H' by 
the conserved variables Q. By combining Eqs. (12) and (14), 
the Euler equations are discretized as 

£7“ I + £ T " V A * ,AQ i + £ t ' 1a " a « AQ m 

2 At m-1 

2 + <|> Q~ -Q" v Q n - Q"- 1 
2 At 2 At 

- X T 1 1 H* (q” ) + H” (q + ) ]* A, (15) 

where I is the identity matrix. Direct solution of the system of 
simultaneous equations which results from application of Eq. 
(15) for all tetrahedra in the mesh, requires the inversion of 
a large matrix with large bandwidth which is computationally 
expensive. Instead, a Gauss-Seidel relaxation approach is used 
to solve the equations whereby the summation involving A Q m 
is moved to the right hand side of Eq. (15). The terms in this 
summation are then evaluated for a given time step using the 
most recently computed values for the aQ's. The solution 
procedure then involves only the inversion of a 5x5 matrix 
(represented by the terms in square brackets on the left hand 
side of Eq. (15)) for each tetrahedra in the mesh. The 
procedure is implemented by first ordering the elements that 
make up the unstructured mesh from upstream to downstream, 
and the solution is obtained by sweeping two times through the 
mesh as dictated by stability considerations. The first sweep 
is performed in the direction from upstream to downstream 
and the second sweep is from downstream to upstream. For 
purely supersonic flows the second sweep is unnecessary. 


The implicit algorithm is formulated by first 
approximating the time derivative in the Euler equations by 

d Q 2 + <j> A Q 2 + <t> Q* - Q n Q n - Q 0 " 1 

at “ 2 At + 2 At 2 At (12) 

where AQ = Q n + 1 -Q* and where the parameter $ controls 
the temporal order of accuracy. For example, the scheme is 
first-order-accurate in time if - 0 and the scheme is 
second-order-accurate in time if $ - 1. For an implicit 
temporal discretization, the flux H must be treated at time 
level (n+1) which is accomplished by linearizing according to 


3 Q 


Q-Q 


AQ 


(13) 


where 3H/9Q is the flux Jacobian A. Also, Eqs. (12) and (13) 

involve q\ the vector of flow variables at an iterate level O 
which is normally taken to be time level (n). For unsteady 
applications, however, subiterations may be performed to 

drive Q to Q n + 1 and thus minimize linearization and 
relaxation errors. 


With the flux-split spatial discretization, the forward and 
backward fluxes are linearized for a given tetrahedron j as 


Boundary Conditions 


To impose the flow tangency boundary condition along the 
surface of the vehicle, the flow variables are set within 
dummy cells that are effectively inside the geometry being 
considered. The velocity components within a dummy cell, (u, 
v, w)d, are determined from the values in the cell j adjacent 
to the surface, (u, v, w)j. This is accomplished by first 
rotating the components into a coordinate system that has a 
coordinate direction normal to the boundary face. The sign of 
the velocity component in this direction is changed (hence 
imposing no flow through the face) and the three velocity 
components are then rotated back into the original x, y, z 
coordinate system. After considerable algebra this yields 


u 


1-2n x -2n x n y -2n x n z 


u ' 

v 

' = 

-2n x rty 1-2n y -2n y n z 


V 

w 

d 

-2n x n z -2n y n z 1-2nf 


w 


(16) 


where n x , n yt and n z are the x, y, and z components of the unit 
vector that is normal to the boundary face. Also, pressure and 
density within the dummy cell are set equal to the values in the 
cell adjacent to the surface. 
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After application of the upwind-biased interpolation 
formula to determine q* and q+ at each face, the velocity 
components are corrected to give a "strong" implementation of 
the surface boundary condition according to 

“corrected = u-n x (un x + vr^ + wn z ) 

Corrected = v - n y ( un x + vn y + wn 2 ) (17) 

Wcorrected = w - n 2 ( un x + vn y + wn 2 ) 

In the far field a characteristic analysis based on Riemann 
invariants is used to determine the values of the flow variables 
on the outer boundary of the grid. This analysis correctly 
accounts for wave propagation in the far field which is 
important for rapid convergence to steady-state and serves as 
a "nonreflecting" boundary condition for unsteady applications. 


BesuHs and Discussion 

To assess the accuracy and efficiency of the unstructured- 
grid upwind-Euler algorithm and to demonstrate several of the 
basic features of the scheme, calculations were first 
performed in two dimensions for the NACA 0012 airfoil. 
These results were obtained using the unstructured grid shown 
in Fig. 2. The grid has 3300 nodes, 6466 triangles, and 
extends 20 chordlengths from the airfoil with a circular outer 
boundary. Also there are 110 points that lie on the airfoil 
surface. Steady-state calculations were performed for the 
airfoil at a freestream Mach number of M_« 0.8 and an angle 
of attack of a » 1.25°. The results were obtained using both 
the implicit relaxation time-marching scheme and the explicit 
four-stage Runge-Kutta time-marching scheme. The explicit 
time-marching results were obtained using a CFL number of 
2.5 (since the CFL limit is approximately 2.8) and the 
implicit time-marching results were obtained using a CFL 
number of 100,000. Such a large value was used for the 
implicit results since the relaxation scheme has maximum 
damping and hence fastest convergence for very large time 
steps. This is in contrast with implicit approximate 
factorization schemes which have maximum damping for CFL 
numbers on the order of 10. 

A comparison of the convergence histories between explicit 
and implicit time-marching is shown in Fig. 3(a). The 
"error" in the solution was taken to be the L 2 ~norm of the 
density residual. As shown in Fig. 3(a), the explicit solution 
is very slow to converge. This solution takes approximately 
10,000 time steps to become converged to engineering 



accuracy, which is taken to be a four order of magnitude 
reduction in solution error. In contrast, the implicit solution 
is converged to four orders of magnitude in only approximately 
500 steps and is converged to machine zero in less than 2000 
steps. The implicit solution costs approximately 75% more 
per time step than the explicit solution because of the 
increased number of operations required to evaluate the flux 
jacobians. This increase in CPU time is far out-weighed by 
the faster convergence to steady state in that a converged 
solution is obtained with the implicit relaxation scheme with 
an order of magnitude less CPU time than the explicit scheme. 
The resulting steady pressure distribution is shown in Fig. 
3(b). For this case there is a relatively strong shock wave on 
the upper surface of the airfoil near 62% chord and a 
relatively weak shock wave on the lower surface near 
30% chord. The pressure distribution indicates that there 
is only one grid point within the shock structure, on either the 
upper or lower surface of the airfoil, due to the sharp shock 
capturing ability of flux-vector splitting. 


log _4 
(error) 



explicit 

time-marching 


implicit 

time-marching 


i L 


5000 


Iteration 


(a) convergence histories. 



(b) steady pressure distribution. 


Fig. 2 Partial view of unstructured grid of triangles about the 
NACA 0012 airfoil. 
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Fig. 3 Comparison of steady-state results for the NACA 0012 
airfoil at M w - 0.8 and a - 1.25°. 


To assess the accuracy of the three-dimensional, 
unstructured-grid, upwind-Euler algorithm, calculations 
were performed for the ONERA M6 wing. 19 The M6 wing has a 
leading edge sweep angle of 30°, an aspect ratio of 3.8, and a 
taper ratio of 0.562. The airfoil section of the wing is the 
ONERA "D" airfoil which is a 10% maximum thickness-to- 
chord ratio conventional section. The results were obtained 
using a grid which has 154,314 nodes and 869,056 
tetrahedra. The surface triangulation for the upper surface of 
the wing is shown in Fig. 4 and a partial view of the mesh in 
the symmetry plane is shown in Fig. 5. The mesh in the 
symmetry plane reveals how coarse the mesh actually is off of 
the surface of the wing. The results were obtained using the 
explicit time-marching scheme since it requires half of the 
memory of the implicit scheme. The code was run for 
6000 time steps at a CFL number of 5.0. The solution 
required approximately 150 hours of CPU time and 125 
million words of memory on the Cray-2 computer at the 
Numerical Aerodynamic Simulation facility located at NASA 
Ames Research Center. 

Figure 6 shows surface pressure coefficient comparisons 
with the experimental data at five span stations including 
q ■ 0.2, 0.44, 0.65, 0.9, and 0.95. In these plots the Euler 
results are given by the solid curves where plus signs have 
been included to indicate the actual grid point values which are 



Fig. 4 Upper surface grid for the ONERA M6 wing. 



Fig. 5 Partial view of symmetry plane mesh for the ONERA 
M6 wing. 


connected with straight line segments. The experimental data 
is represented by the circles. For q « 0.2 shown in Fig. 6(a), 
there are two shock waves along the chord. The forward shock 
wave is well predicted including the suction peak. The second 
shock wave is predicted slightly downstream of the 
experimental shock location which is typical of inviscid 
methods for this case. Also, the lower surface pressure 
coefficients agree well with the data. At q « 0.44 shown in Fig. 
6 (b), the shock locations have begun to coalesce. The leading 
edge suction peak is well predicted and both shock waves are 
captured sharply. At q - 0.65 shown in Fig. 6(c), the forward 
shock wave is near 20% chord and the second shock wave is 
near midchord. All of the pressure levels are well predicted 
and both shocks are captured sharply with only one grid point 
within the shock structure. There are also no overshoots or 
undershoots near the shocks due to the flux limiting. 
Furthermore, the lower surface pressure coefficients are 
predicted accurately. At q = 0.9 shown in Fig. 6(d), the two 
shocks have merged to form a single, relatively strong, shock 
wave near 25% chord. Here the shock is very sharply 
captured and the calculated pressures again agree well with the 
experimental data. Finally at q = 0.95 shown in Fig. 6(e), the 
shock wave is slightly stronger than the previous span station. 
Here, the calculated shock also has only one interior point. 
For this case though, the pressure level upstream of the shock 
is slightly underpredicted due to the coarseness of the mesh. 

Figure 7 shows pressure contour lines on the surface of 
the wing plotted using an increment of ap - 0.02. Pressure 
contours on the upper surface are shown in Fig. 7(a); 
Pressure contours on the lower surface are shown in Fig. 
7(b). The upper surface contours (Fig. 7(a)) clearly show 
the lambda-type shock wave pattern formed by the two inboard 
shock waves which merge together near 80% semispan to form 
the single strong shock wave in the outboard region of the 
wing. Also, the forward shock wave appears to be more 
sharply captured than the second shock wave. This, however, 
is due to the mesh being more dense in this region than in the 
midchord region. The lower surface contours (Fig. 7(b)) 
indicate that there is very little spanwise variation in 
pressure. 


Concluding Remarks 

Improved algorithms for the solution of the three- 
dimensional time-dependent Euler equations are presented for 
aerodynamic analysis involving unstructured dynamic meshes. 
The improvements have been developed recently to the spatial 
and temporal discretizations used by unstructured grid flow 
solvers. The improved spatial discretization involves a flux- 
split approach which is naturally dissipative and captures 
shock waves sharply with at most one grid point within the 
shock structure. The temporal discretization involves either 
an explicit time-integration scheme using a multi-stage 
Runge-Kutta procedure or an implicit time-integration 
scheme using a Gauss-Seidel relaxation procedure which is 
computationally efficient for either steady or unsteady flow 
problems. With the implicit Gauss-Seidel procedure, very 
large time steps may be used for rapid convergence to steady 
state, and the step size for unsteady cases may be selected for 
temporal accuracy rather than for numerical stability. 

Steady flow results are presented for both the NACA 0012 
airfoil and the ONERA M6 wing to demonstrate applications of 
the new Euler solver. The steady results showed that rapid 
convergence to steady state could be achieved with the implicit 
time-marching in comparison with results obtained using 
explicit time-marching. A factor of ten reduction in 
computational cost was obtained for the case that was 
presented. Steady flow results were also presented for the 
ONERA M6 wing to determine the accuracy of the three- 
dimensional capability. Only explicit time-marching was 
used because of the large amount of computer memory that was 
required. The computed surface pressure coefficients showed 
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that the shock waves were captured sharply with only one grid 
point within the shock structure due to the flux-splitting and 
there were no overshoots or undershoots near the shocks due to 
the flux-limiting. The pressure coefficients resembled those 
that are known to be produced by structured-grid Euler 
methods and they also agreed well with the experimental data. 
Future work will focus on validating the implicit time- 
marching procedures in 3D and extending the methods to solve 
the Navier-Stokes equations. 
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